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(-H ■ Abstract. We consider a dilute two-component atomic fcrmion gas with unequal 

populations in a harmonic trap potential using the mean field theory and the local 
density approximation. We show that the system is phase separated into concentric 

C^ . shells with the superfluid in the core surrounded by the normal fermion gas in both 

the weak-coupling BCS side and near the Feshbach resonance. In the strong-coupling 
,.^ BEC side, the composite bosons and left-over fermions can be mixed. We calculate the 

pH , cloud radii and compare axial density profiles systemically for the BCS, near resonance 

O ' and BEC regimes 
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1. Introduction: 

The original BCS state for superconductors considers pairing between two species of 
fermions with equal populations. For a long time, theorists studied the fermion system 
with unequal species, or mismatched Fermi surfaces, and proposed this system may have 
different ground state PP, in particular the so called Fulde-Ferrell-Larkin-Ovchinnikov 
(FFLO) phase. Experimentally, however, such superfluid states remain unclear because 
of the difficulty in preparing the magnetized superconductors. 

Experiments with ultra-cold atoms have opened a new era to study this fermion 
system with unequal populations. Through the Feshbach resonance [2j, the effective 
interaction between atoms can be varied over a wide range such that the ground state 
can be turned from a weak-couphng BCS superfluid to a strong-coupling Bose-Einstein 
condensation (BEC) regime. In the homogeneous system, theoretical studies IH E] 
of the unequal fermion species show that the phase transition must occur when the 
resonance is crossed, in contrast to the equal population case where a smooth crossover 
takes place jHl Ej. Breached pair phase |H|, phase separated states are also proposed 
jHl CD] in this system. 

Two recent experiments fUE] studied the trapped ^Li atoms with imbalanced spin 
populations and obtained the density profiles for various population differences. Both 
groups found the system contains a superfluid core surrounded by a normal fermions 
and provide evidence for phase separation near the crossover. 

In this paper, we study this imbalanced fermion system by mean field approximation 
and evaluate the density profiles for various coupling strengths from weak-coupling 
BCS superfluid to strong coupling BEC regime. In particular, we calculate the axial 
density profiles, superfluid and minority cloud radii, and distinguish between the phase 
separation and Bose-Fermi mixture regimes. This paper is organized as follows: In Sec. 
II we briefly review the mean-field approximation for the dilute fermion atoms with 
unequal populations. In Sec. Ill we present our results for various polarizations from 
weaking-coupling BCS superfluid to strong-couphng BEC side. We show that the axial 
density profiles are constant within the superfluid core and decrease beyond the phase 
boundary for the phase separations but are smoothly decreasing functions for entire 
trap for the mixtures. Finally, we conclude with a briefly summary in Sec. IV. 

While this work was in progress, several theoretical papers have also studied the 
same problem under similar approximations |E1 ^M ^1 E] o^' going beyond pTf IH?] . 
Basically, ref [ISl El 113 HEl UZj also conclude the system is phase separated into 
concentric shells with superfluid in the center and surrounding by leftover fermions 
near the resonance. The strong-coupling BEC limit has also been studied [TH | IT ^ IT^. 
In this case, the composite bosons and unpaired fermions can mix. As the population 
difference increases the unpaired fermions can even penetrate into the superfluid core. 
Our paper provides a more systematic study of the entire BCS-BEC regimes for all 
polarizations. 
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2. Formalism: 

Restricting ourselves to wide Feshbach resonance, the two-component fermion system 
can be described by an effective one-channel Hamiltonian 

k,cr k,k',q 

where ^(j(k) = /i^A;^/2m—yUo- and the index a runs over the two spin components. Within 
the BCS mean field approximation at zero temperature, the excitation spectrum in a 
homogeneous system for each spin is (see e.g. PUI for details) 



where ^o-(k) = h?k'^/2m — fi„ are the quasi-particle excitation energies for normal 
fermions, and — t— i- Fo^' an inhomogeneous system, e.g. the system in a harmonic 
trap, a finite system sized effect should be included [1^]. However, the system can 
be treated as homogeneous locally if the number of particles are sufficiently large. A 
local density approximation, or Thomas- Fermi approximation (TFA), is applied and the 
chemical potential for spin a is replaced by 

f^a{r) = /x° - -mu'^r^ , (3) 

with uj the isotropic trap frequency and r the distance from the trap center. We shall 
show results explicitly only for the isotropic trap. In the local density approximation, 
the densities na{r) depends only on the local chemical potentials /io-(r). Hence the 
density profile of an anisotropic trap can be related to an isotropic one by rescaling the 
spatial coordinates appropriately. We then introduce the average chemical potential 

(4) 





/i(r) ^ 


= ^[f^lir) + 


f^iir) = 


/io - 


X 9 9 

-moj r 


and the difference h 


= P^{r)- 


^i(0]/2 = 


= K- 


/.?)/2. 1 


(J2)) becomes 


^T.i(k, 












r) = ^e(k 


,r)2 + A2(r)T/i 





(5) 

where ^(k, r) = fi}y? j^m — \i{r\ We take spin up to be the majority species so that h 
and E^ are positive always. Then the density profiles in a harmonic trap are 



TT-s (r) = n|(r) + ri|(r) 



(2vr; 
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(6) 



7^/(^t) . (7) 

and the total number of particles A^ = J d^rns{r). Here / is the Fermi function. The 
polarization of the system is defined as 
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Now the pairing field A depends on position also. In the local density 
approximation, it obeys an equation similar to the homogeneous case |^I2UJ: 

■l-f(E^)-f(E^) 



Tfi ./N A / \ I Cl rh 

-A(r) = A r) 



Aira J {271 



E-f + E^ K^k^ 



(9) 



For a given scattering length a, we solve equations (0), © and self-consistently 
for fixed total number of particles A^ and polarization P. The solutions to the "gap 
equation" Q may not be unique. The physical solution is determined by the condition 
of minimum free energy among the multiple solutions. We describe the detail procedures 
in [Appendix A 



3. Results and Discussions 

In this section, we investigate the density profiles for various polarizations and coupling 
strengths from positive detuning BCS superfluid to negative detuning BEC side. With 
the aid of density profiles, we evaluate the radii of the superfluid phase boundaries for 
various cases and compare to the current experimental results. We close this section 
with a discussion of axial density profiles. Phase separation versus Bose- Fermi mixture 
can be clarified through the axial density profiles of the population difference. 

In Fig. we plot the radial density profiles for three different coupling strengths 
1/kpa = —0.61, 0.03, and 2.44 (for different columns) with polarization P = 0.2,0.5, 
and 0.9 (for different rows). The total number of particles are fixed to 2 x 10^. Here 
kp is the Fermi wavevector at the trap center for an ideal symmetric Fermi gas with 
the same total number of particles. In all these plots, the system shows a superfluid 
cloud surrounded by a normal Fermi gas except the case in Fig. ^c) where the system 
is completely in the normal state with the polarization P = 0.9. It is consistent with 
the experimental observation |TT] that the superfluid is destroyed by a sufficiently large 
population difference. We remark here however that the critical population difference 
Pc for the destruction of superfluidity obtained here is much larger than that found 
in the experiment ^T]. For l/kpa = —0.61 here, Pc > 0.6 theoretically whereas an 
extrapolation of the data of [TTj gives Pc < 0.3 (see further discussions below). For less 
polarized system on the BCS side [Fig. (TJ^a) and (b)], it shows a clear phase separation 
between the superfluid and a normal Fermi gas. Note that within the superfluid cloud, 
the population difference is zero and the system is just like the typically unpolarized 
superfluid. Outside the superfluid cloud, both components of the fermions exist in the 
normal state which indicates a fraction of the fermions which are not paired-up even at 
the zero temperature. The density profile of the population difference nd{r) peaks at 
the superfluid phase boundary and decreases gradually toward the edge of the trap. Its 
value is equal to the density profile of the majority when r > r^, the radii of the minority 
cloud. At the superfiuid phase boundary, both the majority and minority density profiles 
exhibit discontinuous which has been observed also by the others |T! H IT ^ ITKj. 

Near resonance [Fig. md)-(f)], similar phase separated states are observed. 
However, most of the minority are paired up in this regime such that the density profiles 
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contains mainly the excess fermions outside the superfluid cloud at all P's. This is 
consistent with the homogeneous normal phase boundary extends near resonance to 
large population difference |lj. Near resonance, the superfluid core survives at P = 0.9 
in our calculations whereas experimentally [llj it vanishes already at P ~ 0.7. 

On the BEG side [Fig. ^g)-(i)], all of the minority are paired up and the excess 
fermions can penetrate into the superfluid core. The system contains a superfluid core 
for any P (if sufficiently deep in the BEG regime, see Fig. |3] below). The system then 
contains three different phases: the purely superfluid, the mixture phase, and the normal 
fermions from the trap center to the edge of the trap. The mixture phase extends toward 
the trap center as the polarization increases. In Fig. ^i) (P = 0.9), the system is highly 
polarized and the excess fermions extend deeply into the trap center. In Appendix B[ 
we give an analytic discussion of the density proflles in the BEG limit. 

The radius r^ of the superfluid core is one of the most interesting quantities in 
current studies of the imbalanced fermion system ^21 El El El EI- From Fig. ^ the 
density proflles of the population difference nd{r) have maxima at the phase boundaries 
between the superfluid and the normal fermions. r^ is thus also the peak position of 
nd{r). In Fig. we plot r^ as function of P for three different coupling strengths, r^ 
behaves quite differently above and below the resonance. On the BGS side the superfluid 
is eliminated when the polarization reaches around 0.65 for l/lkpa) = —0.61 and the 
system becomes completely normal with mismatch Fermi surfaces beyond this critical 
polarization. For large coupling strengths, r^ is flnite except when P is exactly 1, since 
the superfluid is stable for any flnite (7^ 1) polarization in the homogeneous case |3]. 
Except for small (P < 0.1) or large (P > 0.9) polarizations, the sizes of the superfluid 
clouds have maxima near the Feshbach resonance at flxed polarization. 

We also plot the radii (r^) of the minority (spin down fermions) cloud in Fig. El 
Below the Feshbach resonance, this radius is the same as the radius of the superfluid 
core because all the minority of fermions are paired up. However, this won't be true 
above the Feshbach resonance where part of minority are not paired up in this regime. 
Unlike r^, vi decreases monotonically as the coupling strength increases. These two radii 
become identical when the minority are paired up completely when the system reaches 
the BEG regime. 

Due to the experimental constraints, one can not measure the radial density profiles 
directly. Instead the axial density profiles are reported in jT^j. In Fig. 01 we plot the 
normalized axial density profiles na{z) [= J dxdynd{r}] of the population difference for 
different coupling strengths at three fixed polarization P. For the cases with phase 
separations [Figs, [ifa), (b), and (d)-(f)] , the corresponding na{z) are constants for 
z < Tg and have a kink at the phase boundary (for z = Vg). This feature results from 
the population difference nd{r) being zero inside the superfluid cloud (see [Appendix G 



in detail) such that na{z) remains the same value as at the phase boundary 2; = r^. 
For a system with a mixed phase region [Figs, mg)-(i)], na{z) increases smoothly 
toward the trap center even as z < Vg. It reaches a constant value within the cloud 
containing superfluid only [ e.g., the cases with solid lines in Figs. Ufa) and (b)]. At 
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large polarization [solid line in Fig. EI^c)], the excess fermions mix with superfluid entirely 
for 1/kpa = 2.44 such that na{z) increases monotonically toward the trap center. The 
completely different features for the axial density proflles of the population difference 
inside the superfluid cloud can help us to clarify whether the system is phase separation 
or mixture. 

4. Conclusion 

We have investigated the radial density proflles of the two-component fermion system 
with unequal spin-populations under Feshbach resonance. The system shows a 
superfluid cloud in the trap center surrounding by normal fermions. In the weak- 
coupling BCS side, the superfluid is destroyed completely at the polarization P > 0.65 
for l/lkpa) = —0.61. Near the Feshbach resonance, almost all the minority are paired 
up and the system is phase separated into superfluid and the normal fermions. In the 
strong-coupling BEC side, the excess fermions can mix with the superfluid and the 
system contains three different phases, the purely superfluid cloud, the mixture phase, 
and the normal fermions from the trap center to the edge of the trap. 

In particular, we emphasize the difference in the axial density difference profiles 
between the phase separation and Bose-Fermi mixture regimes. The former shows a 
constant for z < Vg and has a kink at the phase boundary but the later are smoothly 
increasing toward the trap center. However, Ref. [T2j reported positive slopes of the 
axial density different profiles inside the superfiuid cloud that indicates the ^^(r) need to 
be negative somewhere within the superfiuid cloud. We do not obtain this phenomenon 
within current local density approach. 
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Appendix A. Free energy 

To determine the minimum free energy state when there are multiple solutions, we need 
an expression for the free energy. We obtain this as follows. First consider a system 
with fixed volume V and particle numbers A'^^-. With the Hamiltonian in Eq. (0), it is 
straight-forward to show that f2I| the energy E{N(j,g) of this system obeys 



dg g 



Asymmetric Fermi superfluid in a harmonic trap 7 

We need to eliminate g in favor of the scattering length a, the physical parameter of the 
system. These two variables are related by the expression 
m 11 ^^ 1 

k 



We thus get 



"4^1=4-1. (A.3) 



and therefore 

dE , T^ "^ I A / m2 /A.N 

Writing this derivative to be £' , we then get the thermodynamic relation 

dE = y" fi^dN^ + £'d(-) . (A.5) 

^-^ a 

a 

The free energy Q{fia, a) = E — ji-^N^ — jiiN^ then obeys 

rff] = - V N„dii„ + £'di-) . (A.6) 

■^^ a 

(7 

We thus conclude that, for volume V and chemical potentials /Xo-, 

Though this expression is already sufficient to determine and thus compare the free 
energies of the different solutions for given scattering length a, we can convert it to an 
even more convenient form. To do this, let us write x = |Ap and y = 1/a. We get, up 
to an overall multiplying factor 

— = -X . A.8 

dy 

Thus, when the solutions are plotted in the form of Fig. the free energy f] of a state 
can be related to the free energy f2o of another state along the same curve by, again up 
to an overall multiplicative constant, 

fi = fio - \ xdy . (A.9) 

Since / xdy is the area between the curve and the y-axis, the states corresponding to 
the minimum free energy at a given scattering length a (hence y) can be found via the 
same procedure as the usual Maxwell construction. 
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Appendix B. BEC limit 

Here we try to understand the behaviour of the density profile in the BEC {1/kpa ^ 1) 
limit, for the special case of small number of excess fermions [e.g. Fig. ^g)]. For a 
bulk system, it is straight forward to perform a low density expansion of the mean-field 
equations and obtain an expansion of the chemical potentials /i/ = /i| and fib = fJ'^+f^i as 
a series in the densities rif = n^ and rif, = n^ (see |22] for the details of this calculation). 
In the BEC limit, n^ can be interpreted as the density of the excess unpaired fermions 
and rif, as the density of the bosons which represent the bound fermion pairs, fif and 
jj^b can be interpreted as the corresponding chemical potentials. Explicitly we have an 
expansion of the form 

f^f = ^^T + 9bfnb , (B.l) 

/ib = - Cfc + gtbnb + gbfUf , (B.2) 

where we have dropped the terms higher order in the densities. These terms are smaller 
in the limit Ufa^ <^ 1 and Uba^ -C 1. We obtain |221 A = (67r^)^/^/2m, e^ = h? /ma^, 
gbb = 4:7ih'^a/m, gbf = 87ih?a/m. The value of A is as expected for a free Fermi gas and 
eb is the binding energy of the fermion pair. The values of gbb and gbf here, obtained 
from the mean-field equations, differ from the correct values known. Exact three-body 
calculation J2SI gives g^}'^'^^ = S.Gnh'^a/m, while a recent four-body calculation |21] gives 
5f™^* = 1.27rh'^a/m. Alternatively, defining the scattering lengths in the usual manner, 
our low density expansion yields the effective values abf = 8a/3 and abb = 2a, whereas 
the exact values should be alY'^^ = l-2a and a^^'^*^* = 0.6a. The difference between the 
mean-field and exact results is due of course to the approximate nature of the mean-field 
theory. We note however, here that (a) abj and abb are both positive and of order a, 
hence in the BEC limit we have necessarily n^a^ -C 1, thus the fermions and the bosons 
can mix J2S1; (b) gbf/dbb > 1/2 for both mean- field and exact values. We shall explain 
the importance of this second relation below. 
In the trap, Eqs. ()B.1|) and ()B.2|) becomes 

fif = Aridirf' + gbMr) + Vir^ , (B.3) 

/ife = - £6 + gbbnb{r) + gbfnf{r) + 2V{r) . (B.4) 

Here V{r) is the trap potential. We shall only consider the case where the potential 
increases from the center of the trap. To appreciate the implications of the relation (b), 
consider the limit of a small number of fermions. The density profile of the bosons can 
then be determined by first ignoring the fermion term in Eq. (jB.3|l . and we obtain the 
boson density profile 

nb{^ = il^b + eb- 2Vir))/gbb , (B.5) 

if the R.H.S. is larger than zero, and Ub = otherwise. For ease of reference, we shall 
call these two regions "inside" and "outside" below. 
Eq. ()B.3|) can be rewritten in the form 

l^f = Anrf(f)2/3 + Ves{f) , (B.6) 
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where Kflf(^) = V{r) + ghf^tjif) is the effective potential for the fermions. We obtain, 
for the inside region, KflF(^) = V{f) + {ghf / gbb){.l^b + ^b — 2V"(r)) whereas for the 
outside Vcsif) = V{r). In the inside region, the position dependence is therefore 
— {2ghf/ghh— l)V{r). From this, it is clear that ii gbf > gbbf^, then the effective potential 
actually is larger near r = than near the edge of the boson cloud. It decreases from 
the center of the trap till it reaches the edge of the boson cloud, then it increases again 
due to the spatial dependence of V{r). Therefore, we conclude that for gbf/gbb > 1/2, 
the excess fermions lie near the edge of the boson cloud, at least for a small number of 
Fermions [2S1E]- The excess fermions lie near the edge of the cloud [Fig. ^g)]. Thus 
a peak in nd{r) does not indicate phase separation. 

Appendix C. Axial Density 

We here discuss some useful expressions governing the axial density within the local 
density approximation. Our results here are, to some extent, a further development of 
those obtained in [Hj. 

In local density approximation, any quantity, say the density difference n^lr), is a 
function entirely of the local chemical potential fi{f) (note that h is a constant). Hence, 
we can write 

nd{r} = g{f^{r}) = g{f^{r}, h) (C.l) 

for some function g. We here note that any density must be identically zero for 
sufficiently large and negative chemical potential, and thus g{C,) vanishes exactly for 
sufficiently large and negative argument C,. It will be convenient to define another 
function G{Q via 

G(C) = r dCgiC) . (C.2) 



Consider now for example the radially integrated density (referred afterwards as 
axial density) difference: 

na{z) = / dxdyridir) . (C.3) 

The integral can be written as, for a trap that is cylindrically symmetric with respect 
to z, 71 Jg°° dp^g{no — ^a^z"^ — |a||P^) where p^ = x^ + y'^. This integral can be expressed 
in terms of G and thus 

27r 1 

na{z) = — G(/io - -c^zz^) . (C.4) 

a\\ 2 

This relation can be used to obtain directly the axial density in terms of the 
function g, without first calculating the density profile in trap and then performing the 
integration. Moreover, it can also be used to deduce the (unintegrated) density profile 
once the axial density is given, for example from experiment. To see this, we differentiate 
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Eq. ()C.4|) with respect to z, and notice that G"(C) = 5'(C) evaluated at ( = fiQ — \oizZ^ 
is directly related to the corresponding density at the coordinate (0, 0, z). Thus we have 

MO.O..) = --^i^. (C.5) 

2'naz z oz 

Therefore the actual density for a point on the z axis can be obtained from the axial 
density. Under the local density approximation, the density at any given point in the 
trap is a function of the local chemical potential only. Hence we can obtain the actual 
density at any given point in space provided the axial density is given. 

The relation Eq. ()C5|) also shows that, since n^ is non-negative, the axial density 
profile is strictly decreasing (increasing) with z for z > (<) 0, a result pointed out in 
reference jHj. In a region where the density vanishes, then the axial density should be 
a constant, a result recognized in references jHj and |15j . 
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Figure 1. (color online). Radial density profiles for n|(r) (solid lines), ni{r) (dashed 
lines), and nd{r) (solid circles). The dotted lines arc the normal state density profile 
with the same density as in the normal region. The coupling strengths (each column) 
l/ikpa) = -0.61 (left, BCS-side), 0.03 (middle, near resonance), 2.44 (right, BEC-side) 
and the polarizations (each row) P = 0.2, 0.5, and 0.9. The total number of fermions 
are fixed to A^ = 2.0 x 10^. Rtf is the size of the Thomas- Fermi radius of the normal 
spin-up cloud. 
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Figure 2. (color online). Radii of the superfluid core Vg versus polarization P for 
normalized coupling strengths in the legend. 
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Figure 3. (color online). Radii of the minority species r^ versus polarization P for 
normalized coupling strengths in the legend. 
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Figure 4. (color online). Axial population density difference na{z) versus z/Rtf- 
Rtf is defined in the caption of Fig 1. na{z) is normalized to 1 at z = r^. 
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Figure 5. Solution to the gap equation plotted as 1/avs |Ap (arbitrary units) for tlie 
case witlr multiple solutions (dashed lines). The minimum energy solution is shown as 
the thick black line. The areas labeled by A and B are equal. 



